
function [psdx_array, wave_psd_avg] = waveno_only_spectra_shear_ybp(input_u,input_v,z_ip,NFFT_waveno_shear)


% Author: Ritabrata Thakur, 1 October 2021
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

win_space = hann(NFFT_waveno_shear);

%Calculate the scaling factor.
S2 = sum(win_space.^2)/NFFT_waveno_shear;

clear data_u data_v

if (length(size(input_u))>2)
  data_u = squeeze(input_u(1,1,:,:));
  data_v = squeeze(input_v(1,1,:,:));
else
  data_u = squeeze(input_u(:,:));
  data_v = squeeze(input_v(:,:));
end


no_psds=1;
    for j=1:length(data_u(1,:))
        z_ip_portion  = z_ip;
        data_portion_u = data_u(:,j);
        data_portion_u = data_portion_u - mean(data_portion_u);%removing the mean         
        %from individual portion of time series
        data_portion_u = detrend(squeeze(data_portion_u), 'omitnan');%removing any trend from 
        %individual portion of time series
        %data_portion_w = whitening(data_portion', Fs, 'freq', [0,Fs/2]); %for the full range of 0 Hz to the Nyquist frequency
        %data_portion_w = data_portion_w';
        data_portion_v = data_v(:,j);
        data_portion_v = data_portion_v - mean(data_portion_v);%removing the mean 
        %from individual portion of time series
        data_portion_v = detrend(squeeze(data_portion_v), 'omitnan');%removing any trend from     
        %individual portion of time series
        
        
        %% Applying spatial filter after removing mean
        %for i=1:length(data_portion_u(1,:))
        %  data_portion_u(:,i) = data_portion_u(:,i).*win_space';
        %end
        
        data_portion_u = squeeze(data_portion_u).*win_space;
        data_portion_v = squeeze(data_portion_v).*win_space;
        %data_portion_u = detrend(data_portion_u, 'omitnan');
        %for i=1:length(data_portion_v(1,:))
        %  data_portion_v(:,i) = data_portion_v(:,i).*win_space';
        %end
        
             
        %% Making the complex velocity signal and calculating shear
        complex_vel = data_portion_u + 1i*data_portion_v;
        
        % Shear according to Pinkel 2014 (pp 2078).
        for ii=1:length(data_portion_v)-1
            %for kk=1:length(data_portion_v(1,:))
              complex_shear(ii) = (complex_vel(ii)-complex_vel(ii+1))./(mean((z_ip_portion(ii)-z_ip_portion(ii+1)), 'omitnan'));
            %end
        end
        
        
        %% calculate one-D discrete Fourier transform on the complex windowed data
        % The number of FFT points is different because shear has one point less than
        % velocities
        xdft_space = fft(complex_shear, length(complex_shear));
        
        % with NFFT number of points
        psdx = xdft_space.*conj(xdft_space); % This is same as abs(xdft).^2;
      
        % For the following normalisation, refer to
        % https://holometer.fnal.gov/GH_FFT.pdf or holometer.pdf section 9
        % (equation 24). The factor 2 is multiplied in the next step

        psdx = psdx/S2 ; %normalising by the weight of the windows.  
        
        psdx = psdx(1:floor(length(complex_shear)/2+1)); 

        
        %% One-sided spectra for real DFT. 
        psdx(2:end-1) = 2*psdx(2:end-1); %multiplying all frequencies except 0 and N/2 by 2 as 
        %they are overlapped on top of each other. It is defined as two times the two-sided PSD 
        %It does not really matter if you multiply the zero and N/2 with 2 or not, we get the
        %same spectra. Also, we have removed the mean so the component in
        %zero frequency does not matter for the total energy
        
        % saving psd from each data portion in a 3D array which will be
        % summed element wise and then the average psd calculated
        psdx_array(:,:,no_psds) = psdx;
        clear data_portion_u data_portion_v
        no_psds = no_psds+1;
        
    end
    %% Averaging the individual and presumably independent estimates
    % The individual psds are not good estimates of PSDs. They contain too
    % much statistical oscillation and hence they are averaged
    % (https://www.osti.gov/servlets/purl/5688766/) or
    % welch_otis_solomon.pdf. psd_avg is the Welch estimate and matches the
    % estimate by signal.welch estimate in Python
    wave_psd_avg = sum(psdx_array,3)/(no_psds-1);
    
end
